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We study the variation of the spectrum of the Fermi Bubbles with Galactic latitude. Far from the 
Galactic plane (|&| ^ 30°), the observed gamma-ray emission is nearly invariant with latitude, and is 
consistent with arising from inverse Compton scattering of the interstellar radiation field by cosmic- 
ray electrons with an approximately power-law spectrum. The same electrons in the presence of 
microgauss-scale magnetic fields can also generate the the observed microwave "haze" . At lower lat- 
itudes (|6| <J 20°), in contrast, the spectrum of the emission correlated with the Bubbles possesses a 
pronounced spectral feature peaking at ~l-4 GeV (in E^dN/dE) which cannot be generated by any 
realistic spectrum of electrons. Instead, we conclude that a second (non-inverse-Compton) emission 
mechanism must be responsible for the bulk of the low-energy, low-latitude emission. This second 
component is spectrally similar to the excess GeV emission previously reported from the Galactic 
Center (GC), and also appears spatially consistent with a luminosity per volume falling approxi- 
mately as r^^ *, where r is the distance from the GC. Consequently, we argue that the spectral 
feature visible in the low-latitude Bubbles is most likely the extended counterpart of the GC excess, 
now detected out to at least ~2-3 kpc from the GC. The spectrum and angular distribution of the 
signal is broadly consistent with that predicted from ~10 GeV dark matter particles annihilating 
to leptons, or from ~50 GeV dark matter particles annihilating to quarks, following a distribution 
similar to, but slightly steeper than, the canonical Navarro-Frenk- White (NFW) profile. We also 
consider millisecond pulsars as a possible astrophysical explanation for the signal, as observed mil- 
lisecond pulsars possess a spectral cutoff at approximately the required energy. Any such scenario 
would require a large population of unresolved millisecond pulsars extending at least 2-3 kpc from 
the GC. 
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I. INTRODUCTION 

Data from the Fermi Gamma-Ray Space Telescope 
have revealed a pair of large gamma-ray lobes extend- 
ing approximately 50° north and south of the Galactic 
Genter [T]. These lobes, known as the Fermi Bubbles, 
are visible in gamma-rays between ~1-100 GeV and have 
a markedly harder spectrum {dN/dE cx E^^) than the 
gamma-ray emission associated with the Galactic Disk. 

The Bubbles were originally studied as a possible 
gamma-ray counterpart to the WMAP haze a spec- 
trally hard microwave excess in the inner Galaxy most 
clearly visible in WMAF''s 23 and 33 GHz frequency 
bands. The haze was first discovered in 2003 [2, and has 
been studied over the past decade as a possible signature 
of a new hard electron population in the inner Galaxy 
[IHS], producing microwave synchrotron radiation in the 
Galactic magnetic field. Recently, the existence of the 
microwave haze has been confirmed by the Planck exper- 
iment [7], whose data indicate a strong degree of spatial 
coincidence between the microwave haze and the gamma- 
ray Bubbles, further supporting the hypothesis that these 
signals have a common origin. Perhaps the simplest pos- 
sibility is that the gamma rays arise from inverse Comp- 
ton scattering (ICS) by the same hard electron popula- 
tion that produces the haze via synchrotron. 



The question of the origin and nature of the Fermi 
Bubbles has been the subject of much debate. One 
key question is whether these gamma-rays are produced 
by a hadronic [SHTU] or leptonic [TT| - fT5] mechanism, i.e. 
whether they arise from the scattering of energetic pro- 
tons on the gas of the interstellar medium, or from the 
ICS of photons from the interstellar radiation field by en- 
ergetic electrons. An example of a hadronic scenario was 
proposed by Aharonian and Crocker [5], in which the 
Bubbles are billion-year-old reservoirs of energetic pro- 
tons, which were injected as a result of star formation 
in the Galactic Genter and are confined within the Bub- 
bles by magnetic fields. Leptonic scenarios have garnered 
somewhat more attention in the literature, and provide 
a straightforward link with the spatially correlated emis- 
sion observed in the microwave and radio [16j (hadronic 
scenarios will also generate synchrotron emission through 
the electrons produced in charged pion decays, but since 
these electrons will diffuse after being produced, the con- 
nection between their spectrum and spatial distribution 
and that of the gamma-rays is not as straightforward as in 
the leptonic scenario). In such leptonic models, electrons 
are accelerated by a shock or a series of shocks and/or 
by Fermi acceleration in turbulent magnetic fields behind 
the shock front [17]; the shock(s) may be fueled by ac- 
cretion onto the supermassive black hole at the Galactic 
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Center, by starburst activity, or by some other mecha- 
nism. 

In this article, we reexamine the Fermi Bubbles and 
the variation of their spectrum with Galactic latitude. 
Far from the Galactic plane (|&| > 30°), the observed 
gamma-ray spectrum is nearly invariant with latitude 
and fairly flat over the energy range observed by Fermi. 
This spectrum can be well explained by inverse Comp- 
ton scattering of cosmic microwave background (CMB), 
infrared, and starlight photons by a population of GeV- 
TeV electrons with an approximately power-law spec- 
trum {dNe/dEe ^ E~^). Furthermore, we find that this 
same population of cosmic ray electrons leads to syn- 
chrotron emission of the same amplitude as the observed 
microwave haze, if microgauss-scale magnetic fields are 
present in the high-latitude regions of the Bubbles. The 
success of this simple and self-consistent picture provides 
strong support for a leptonic origin of the high-latitude 
emission from the Fermi Bubbles. 

At latitudes closer to the disk, however, a leptonic ori- 
gin of all the emission associated with the Bubbles does 
not appear possible. The gamma-ray spectrum of the 
Fermi Bubbles at latitudes within approximately 20° of 
the Galactic plane possesses a peak at energies of a few 
GeV, and cannot be generated by inverse Compton scat- 
tering of starlight, infrared, or cosmic microwave back- 
ground radiation by any realistic steady-state electron 
population. Furthermore, no realistic spectrum of cosmic 
ray protons is capable of accounting for the gamma-ray 
spectrum observed at these low latitudes. 

Gamma-ray emission with a similar spectrum has been 
previously identified from the region surrounding the 
Galactic Center (GC) Proposed origins for this 

excess include annihilating dark matter [T8H201 [22], a 
population of millisecond pulsars [THl Uni [^2HM] , or cos- 
mic ray interactions with gas [T51 [121 [HI US [5B] . In this 
paper, we show that the non-inverse Compton compo- 
nent of the emission from the Fermi Bubbles identified 
in this study is spectrally and morphologically consistent 
with being the extended counterpart of this GC excess, 
revealing that this emission is not confined to the GC, 
but extends out to at least ^2-3 kiloparsecs from the 
Galactic plane. The morphology of this signal is consis- 
tent with originating from annihilating dark matter dis- 
tributed according to a generalized Navarro-Frenk- White 
(NFW) profile with an inner slope of p cx r~^'^, where 
r is the distance to the GC (the GC excess, in isola- 
tion, favors a power-law slope in the range of 1.2-1.4 
[m [HI [22] ) . The spectral shape of this signal can be ac- 
commodated by dark matter particles of mass ~10 GeV 
annihilating to leptons, or by ~50 GeV particles annihi- 
lating to quarks. In either case, the normalization of the 
observed signal requires an annihilation cross section on 
the order of ctw ~ (6 — 8) x 10"^'' cm'^/s, similar to that 
expected of a thermal relic of the Big Bang. 

The remainder of this paper is structured as follows. 
In Sec. [n] we describe the analysis used to extract the 
spectra from various regions of the Fermi Bubbles. In 



Sec. Ill we demonstrate that the high-latitude regions of 
the Bubbles can be accounted for with an approximately 
power-law spectrum of GeV-TeV electrons, which can si- 
multaneously produce the observed microwave haze as 
synchrotron. In Sec. IV we turn our attention to the 



low-latitude regions of the Fermi Bubbles, and show that 
their spectra cannot be accounted for by inverse Comp- 
ton scattering. Instead, an additional mechanism is re- 
quired, capable of producing a spectrum which peaks 
strongly at energies of a few GeV. In Sec. |V] we com- 
pare this signal with that previously observed from the 
Galactic Center. In Sec. |VI[ we study the residuals re- 
maining when the (best-fit) known backgrounds are sub- 
tracted, and demonstrate how the few-GeV peak appears 
in these residuals. In Sec. |VII[ we show that the emis- 
sion observed at low-latitudes is better described by a 
spherically symmetric, NFW- like morphology than by a 
flat-brightness distribution confined to the regions of the 
Fermi Bubbles. In Sec. |VIIl"] we discuss possible interpre- 
tations of this signal, including annihilating dark matter 
and a population of gamma-ray pulsars. We summarize 
our results and draw conclusions in Sec. |IX[ This paper 
also includes five appendices, which describe various cross 
checks of our results and other supplementary material. 



II. EXTRACTING THE 
LATITUDE-DEPENDENT SPECTRUM OF THE 
GAMMA-RAY BUBBLES 

In our analysis, we employ the publicly available Pass 
7_Version 6 data release from Fermi, with 4.5 years of 
photon data.^ We apply a standard zenith angle cut 
to exclude the Earth limb, rejecting events with zenith 
angles greater than 100°. We also employ the recom- 
mended diffuse analysis cuts on the data quality, nominal 
science configuration and rocking angle: DATA_QUAL==1, 
LAT_C0NFIG==1, ABS (ROCK_ANGLE) <52. Throughout, we 
use the class of events designated ULTRACLEAN, but have 
confirmed that employing the SOURCE or CLEAN classes 
does not alter our conclusions. 

We generate skymaps for 30 log-spaced energy bins 
spanning the range from 0.3 GeV to 300 GeV, binning the 
photons on an NSIDE=256 HEALPix grid, and smooth all 
maps to 2° FWHM (full width at half maximum). As in 
Ref. |T], we use front-converting events only (which have 
better inherent angular resolution) at energies below 1 
GeV. We follow the prescription from Ref. [1] for point 
source subtraction, using the Fermi 1-year source catalog 
and masking the very brightest sources. 

For each energy bin, we fit the skymap as a linear com- 
bination of templates, maximizing the Poisson likelihood. 
The Gaussian errors on the fit coefficients are computed 



^ The dataset we employ may be downloaded from 
http: //heasarc . gsfc .nasa.gov/FTP/fermi/data/lat/weekly. 
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FIG. 1: The regions of the sky considered in our analysis. 
The Fermi Bubbles themselves are broken into five pairs of 
regions according to Galactic latitude (|6| < 10°, 10° — 20°, 
20°-30°, 30°-40°, and 40°-50°). Also shown as dashed lines 
is the inner complement region to the Bubbles, as described 
in Appendix [P] 

from the likelihood by Aln£ = 1/2, and do not take 
into account the systematic error in the event that the 
linear combination of templates is a poor description of 
the data. Further details of the fitting procedure may be 
found in Ref. [1] and in Appendix B of Ref. . We em- 
ploy several different template combinations to test the 
robustness of our results to the foreground model. 

In the Galactic disk, there is a substantial popula- 
tion of unsubtracted point sources, as well as bright 
diffuse emission; consequently, we mask the inner disk. 
Throughout our study, we will show results for masks 
with |6| < 1°, |6| < 2° and |&| < 5°, to test the depen- 
dence of our results on this parameter. 

To determine the spectrum of the Fermi Bubbles as a 
function of latitude, we divide the standard spatial tem- 
plate for the Bubbles (as defined in Ref. pj) into sub- 
regions by (absolute) latitude: |6| < 10°, 10° < |6| < 20°, 
20° < |5| < 30°, 30° < l&l < 40°, and 40° < |&| < 50° 
(see Fig. [T]). We smooth all templates to the scale of 
the maps. We fit separately for the spectrum in each of 
these latitude bands, varying the degree of masking of 
the Galactic Disk, and with a range of template-based 
models for the known backgrounds. All fits include an 
isotropic offset, to absorb residual cosmic-ray contami- 
nation and isotropic diffuse emission. The two primary 
possibilities we consider for the additional templates are: 

• Diffuse model: We take the Fermi diffuse model 
from version P6V11 of the Fermi Science Tools, 
smooth it to match the maps we are using, in- 
terpolate to the appropriate energies, and perform 
the fit using only this template (in addition to the 



latitude-sliced Bubbles templates and the isotropic 
template, which are universal to all subtraction 
methods). This version of the diffuse model has 
been adjusted to fit the data assuming no contri- 
bution from the Bubbles; consequently it may ab- 
sorb some of the Bubbles-correlated emission at the 
cost of oversubtraction in neighboring regions. It 
was also designed primarily to model the emission 
at energies < 50 GeV, and is not recommended for 
use at very low latitudes, |&| < 3°. However, since 
our signal extends to quite high latitudes and the 
energies of greatest interest are at < 50 GeV, these 
latter caveats do not pose severe problems for our 
study. 

• Low- energy template: We employ the Schlegel- 
Finkbeiner-Davis (SFD) dust map [37] as a tem- 
plate for emission from cosmic-ray protons scat- 
tering on the gas (see Refs. [3 Hj for a discus- 
sion). We take the Fermi data at 0.5-1.0 GeV 
(where the Bubbles are less pronounced ^) and 
subtract the SFD dust map to obtain an approxi- 
mate template for emission from inverse Compton 
scattering by cosmic-ray electrons; this is the dom- 
inant contribution to the diffuse background after 
the dust /gas-correlated emission has been removed. 
We then fit the higher-energy data using this tem- 
plate, the SFD dust map, and a flat template for 
the large soft-spectrum structure known as Loop L 
This method avoids the use of complicated mod- 
els and minimizes the use of external maps, but by 
construction cannot probe the Bubbles spectrum 
at energies around or below 1 GeV, and does not 
take into account spectral variation in the various 
emission components with position in the Galaxy. 

Each of these templates has been discussed in greater 
depth in Ref. [1]; we refer the reader to that work for 
further details. We also employ the same normalization 
convention as Ref. [T]; the coefficients of the SFD dust 
map, 0.5-1.0 GeV map, and Fermi diffuse model are mul- 
tiplied by the average value of these templates within the 
entire region defined by the Bubbles (and outside of the 
mask). The other templates are flat (in projected inten- 
sity) within the regions that they are non-zero (over a 
given latitude range within the Bubbles, for example). 

We employ the "diffuse model" fit for our primary re- 
sults, and use the "low-energy template" fit as a check 
for possible systematic errors introduced by our use of 
the Fermi diffuse model. We have also tested the "sim- 
ple disk" model as employed in Ref. where the ICS 
emission is described by a simple geometric diffused disk 
template. As this template proved inadequate for mod- 
elling the data close to the Galactic plane, and the results 
obtained using it were found to depend strongly on the 
degree of masking of the disk, we have relegated discus- 
sion of this model to Appendix [A] 

In Fig. [2] we show the extracted spectrum for each of 
the fitted templates, masking only the region within one 
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FIG. 2: The spectra of the various fit components, including five separate latitude-shced templates for the Bubbles (see Fig.[T|, 
for the two foreground models we employ (see text). The Galactic Disk is masked for |6| < 1° in each case. The left and center 
panels employ the "diffuse model" fit, for the entire sky in the left panel and the southern hemisphere in the center panel. The 
right panel employs the "low-energy template" fit over the entire sky (see text for the details of the fitting procedures). 



degree of the Galactic plane, |&| < 1°. We show results 
found using the "diffuse model" and the "low-energy tem- 
plate" . In the center frame, we show the fit restricting 
to the southern sky (5 < 0), which we might expect to 
be less contaminated by bright features such as Loop I. 
As expected, the error bars are larger in this case due to 
lower statistics, but the results are not otherwise signifi- 
cantly altered. In Fig. |3j we show the spectra extracted 
for the gamma-ray Bubbles, and the dependence on the 
degree of masking of the disk, in each range of Galac- 
tic latitude. For our two preferred template models, the 
results are largely stable to changes in the mask. 

While the gamma-ray spectra extracted using the low- 
energy template appear somewhat different from those 
derived using the diffuse model, this is natural and ex- 
pected, particularly at low energies, since part of the 
emission associated with the Bubbles is included in the 
low-energy template itself. Additionally, at high energies 
and low latitudes the low-energy template fit yields a 
significant amount of emission roughly flat in E^dN/dE, 
which is nearly absent in the diffuse model fit; these is- 
sues are discussed in further detail in Appendix [B| 

We note that the spectrum is almost invariant from 
|5| = 20° — 50°. This suggests that the electrons respon- 
sible for the observed emission in any leptonic scenario 
must either be accelerated in situ or instead travel from 
the inner Galaxy very rapidly, avoiding significant energy 
losses (the distance over which TeV electrons propagate 
via standard diffusion without significant energy losses is 
considerably less than the 5 or more kpc to which this 
angular range corresponds). In contrast, a pronounced 
change in the Bubbles' spectrum is observed at lower lat- 
itudes. In an attempt to quantify the significance of this 
transition, we have compared the quality of the fit found 
using five separate latitude-sliced Bubbles templates to 
that found using only a single Bubbles template. Even 
conservatively limiting our analysis to the cleaner south- 



ern bubble, and masking within 5° of the disk, we find 
that the five-Bubbles-templates model is favored over the 
single Bubbles template at the level of approximately 
16(7. However, it is important to note that this is a, formal 
significance, accounting only for statistical error; there is 
a degree of unavoidable and unaccounted-for systematic 
error in that neither model is a "good fit" , in the sense 
of describing the sky to the level of Poisson noise. 



III. COSMIC RAY ELECTRONS AS THE 
SOURCE OF THE HIGH-LATITUDE 
GAMMA-RAY BUBBLES AND SYNCHROTRON 
HAZE 

Following Blumenthal and Gould [28], we employ the 
full Klein-Nishina formula to compute the spectrum of in- 
verse Compton emission from an arbitrary electron pop- 
ulation. For the problem at hand, we need to consider 
scattering with the CMB as well as with starlight and 
infrared radiation. In our calculations, we adopt the in- 
terstellar radiation model of Ref. . At energies below 
~3 X 10~^ eV, the CMB dominates the energy density, 
while starlight is important at higher energies. 

The gamma-ray spectra observed from various regions 
of the Fermi Bubbles are shown again in Fig. |4] (as 
found using the diffuse model template fit, and mask- 
ing within 1 degree of the disk). To determine whether 
these gamma-rays could be the product of inverse Comp- 
ton scattering, we take an arbitrary (binned) spectrum 
of electrons and compare the resulting inverse Compton 
emission to that shown in Fig. |4| In Fig. [5] we plot the 
electron spectrum which provides the best possible fit to 
the gamma-ray spectrum for each latitude range (and 
error bars around the best fit). The solid line in each 
frame of Fig. |4] denotes the best-fit spectrum of inverse 
Compton photons. At high latitudes, an approximately 
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FIG. 3: The spectrum extracted for the gamma-ray Bubbles in ten-degree latitude bands: in order from the top row, 
40° < |61 < 50°, 30° < |b| < 40°, 20° < |6| < 30°, 10° < |&| < 20°, |6| < 10°. The left and center panels use the "diffuse 
model" template fit (see text); in the center panels, the fit is restricted to 6 < in addition to the masking. The right panels 
use the "low-energy template" approach (see text). The different colors show different choices for the latitude cut to remove 
the Galactic Disk: |&| < 1° (black), jfoj < 2° (blue), |6| < 5° (red). Where the la error bars overlap with zero, we instead plot 
downward-pointing arrows corresponding to the 3(7 upper limits on the emission. 



6 



3 - 



10" 



> 

CD 

o 



3 - 



10" 



3 - 



10" 




0.5 1.0 5.0 10.0 50.0100.0 

E.^ (GeV) 



10 



-6 



6 

o 

> 

0) 

H 



10 



5 - 



-7 _ 



5 - 



10" 




0.5 1.0 5.0 10.0 50.0100.0 

E, (GeV) 



10" 



> 
w 

T3 

\ 

:z 

Xl 



10" 



10" 




0.5 1.0 5.0 10.0 50.0100.0 

E, (GeV) 



10 



-6 



> 
H 

-d 



5 - 



2 4. 



10 



-7 



5 - 



10" 




0.5 1.0 



5.0 10.0 50.0 100.0 

E^ (GeV) 



10 



-6 



> 
xl 



\ 5 - 



10 



-7 



10" 




0.5 1.0 5.0 10.0 50.0 100.0 

E^ (GeV) 



FIG. 4: The gamma-ray spectrum of the Fermi Bubbles, broken into different regions by Galactic latitude (see Fig. [T]). The 
solid lines denote the best-fit spectrum of inverse Compton emission, as calculated from the central values of the electron spectra 
shown in Fig. [5] At high latitudes, the spectrum of the Fermi Bubbles is consistent with originating entirely from the inverse 
Compton scattering of GeV-TeV electrons, while at lower latitudes inverse Compton scattering alone cannot account for the 
observed emission. The dashed line in each frame denotes the spectrum of inverse Compton scattering that would be predicted 
from a spectrum of electrons the same as that required to generate the inverse Compton scattering spectrum observed in the 



highest latitude region (|&| = 40° — 50°), as discussed in Sec. IV 
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FIG. 5: The spectrum of electrons within the volume of the Fermi Bubbles that is best able to fit the observed gamma-ray 
spectrum (see Fig. |4| through inverse Compton scattering alone. At high latitudes, the electron spectrum yields a good fit to 
the observed gamma-ray spectrum, while at low latitudes no spectrum of electrons is able to produce a gamma-ray spectrum 
consistent with observations. 



power-law spectrum of electrons {dNe/dEf, cx E^'^) be- 
tween GeV and TeV energies can produce a spectrum 
of gamma-rays consistent with that observed from the 
Fermi Bubbles. At |6| = 40° - 50°, for example, the 
best-fit electron spectrum shown in Fig. [5] provides an 
excellent fit (x^ = 21.8 over 29-7 degrees-of-freedom) to 



the observed gamma-ray spectrum. 

As noted in the Introduction, the inverse Compton in- 
terpretation of the Fermi Bubbles is supported by the 
observation of spatially correlated emission at both mi- 
crowave j7j and radio 116] wavelengths. For a given 
magnetic field strength (or distribution of magnetic field 
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FIG. 6; The spectrum of synchrotron emission predicted from the best-fit electron spectra as shown in Fig. [S] for various 
values of the magnetic field (see legend). Each frame represents a range of latitudes of the Fermi Bubbles: 40° — 50° away from 
the Galactic plane [top left), 30° — 40° away from the Galactic plane [top right), 20° — 30° away {bottom left), and 10° — 20° 
away {bottom right). Red data points indicate the spectrum of the haze in the southern sky within each given latitude range, as 
derived from Ref. [6]; the average over I is performed for —5° < I < 15°, as in Fig. 5 of Ref. [6]. At somewhat lower frequencies, 
the green triangles and dashed lines connecting them denote the spectral index {not amplitude) for polarized emission between 
2.3 GHz and 23 GHz, as extracted from S-PASS |30] data by Ref. [1^; the averaging is performed over —5° < I < 15° and 
restricted to |&| > 15° to avoid depolarization associated with the Galactic plane. This comparison demonstrates that in the 
presence of microgauss-scale magnetic fields, the electron population responsible for the observed ICS emission can also account 
for the observed synchrotron haze. Also shown for comparison as blue data points is the spectrum of the microwave haze as a 
whole, as given by Ref. [7]. As these results are not binned by latitude, however, they should not be directly compared to the 
predicted synchrotron spectra. 
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strengths), we can calculate the spectrum of synchrotron 
emission that is predicted to result from the spectrum 
of electrons shown in Fig. [5] For simplicity, we take the 
central values for the extracted spectrum and assume an 
isotropic electron population; we do not extrapolate to 
energies above the highest bin. In Fig.[6j we show the re- 
sulting synchrotron spectra for a range of magnetic fields 
from 0.1-100 /xG, assumed to be uniform throughout the 
region of the Bubbles template in question. For compar- 
ison, the magnetic field strength in the local region of 
the Milky Way is thought to be of order a few fxG [3T] . 
One should keep in mind, however, that magnetic fields 
in localized regions and filaments may be much higher 
than the spatially averaged value. 

To compare the predicted synchrotron spectrum to ob- 
servations, we employ the WMAP7 results of Ref. [B], and 
the spectral index of the polarized emission between 2.3 
GHz and 23 GHz, as obtained from S-PASS [3D] data 
and described in Ref. |16j.^ We also show the recent re- 
sults combining data from WMAP and Planck [7|, but 
these results were not binned by latitude so should not be 
directly compared to the predicted synchrotron curves. 

The amplitudes of the gamma-ray Bubbles and mi- 
crowave haze appear consistent with arising from a com- 
mon population of cosmic ray electrons for magnetic 
fields on the order of a few ^G. The spectrum deduced 
from WMAP is somewhat harder than expected from 
the corresponding gamma-ray data, but this may be the 
result of contamination by the CMB (which would be ex- 
pected to produce the appearance of spectral hardening) . 
In contrast, the spectral index for the polarized emission 
found by the authors of Ref. [16] appears consistent with 
expectations from the gamma-rays. If the spectral hard- 
ness observed in WMAP data is confirmed by Planck, 
however, it could suggest a scenario in which most of the 
synchrotron emission is produced in regions with higher 
than average magnetic fields. For example, if magnetic 
fields as high as ^ 30/iG are present in even ^1% of the 
volume of the Bubbles, then this discrepancy could be 
ameliorated. 

We caution that the success of this picture does not 
rule out hadronic scenarios; we focus on leptonic scenar- 
ios in this work because their consistency can be checked 
in a straightforward and model-independent way, since 
the microwaves and gamma-rays originate from the same 
steady-state electron population. In hadronic scenarios, 
in contrast, the gamma-rays provide a probe of the pro- 
ton CR spectrum, but the microwaves probe the electron 
spectrum after diffusion and cooling, and so the consis- 
tency may depend on diffusion parameters in addition 
to the magnetic field. Consequently, we leave a careful 
study of consistency in the hadronic case for future work. 



^ We thank Greg Dobler for providing us with the numerical values 
corresponding to Fig. 5 of Ref. |0 and the authors of Ref. [16] for 
providing the spectral index map corresponding to their Fig. S4. 



Later in this article, we treat the inverse-Compton-like 
emission present in the Bubbles as a background to be 
subtracted; this "background" is an approximately flat 
spectrum in E'^dN/dE with a downturn around 50 GeV, 
and while it is modeled as inverse Compton emission, 
instead treating it as emission from a hadronic cascade 
should not impact our results in any significant way. 



IV. EVIDENCE OF NON-INVERSE COMPTON 
EMISSION AT LOW LATITUDES 

In the previous section, we demonstrated that the spec- 
trum observed from the high latitude (|5| > 30°) regions 
of the Fermi Bubbles can be accounted for by inverse 
Compton scattering of an approximately power-law spec- 
trum of cosmic ray electrons. The same electrons, in 
the presence of /xG-scale magnetic fields, can also ac- 
count for the observed amplitude of the WMAP haze. 
At lower Galactic latitudes, however, we find that in- 
verse Compton scattering alone cannot account for the 
observed gamma-ray emission. In particular, we find that 
the gamma-ray spectrum climbs rapidly between 0.3 and 
1 GeV at low-latitudes (see Fig. |4|, and this rise cannot 
be accounted for by the inverse Compton scattering of 
any physically realistic spectrum of electrons. Quantita- 
tively, for the choice of binning used in Fig. [5| we find 
that an entirely inverse-Compton origin for the gamma- 
rays observed from the \b\ = 1°-10° or |6| = 10° -20° re- 
gions of the Fermi Bubbles yields best-fits of = 100.4 
and 95.7, respectively, each over 29-7 degrees-of-freedom. 
From this, we are forced to conclude that a non-negligible 
fraction of the gamma-ray emission observed from the 
low-latitude regions of the Fermi Bubbles is not the re- 
sult of inverse Compton scattering.^ 

If we assume that the spectral shape of electrons 
present throughout the volume of the Fermi Bubbles 
does not change significantly with latitude, and that 
the gamma-ray emission from the highest latitude region 
(|5| — 40° — 50°) is dominated by the products of in- 
verse Compton scattering, we can calculate the inverse 
Compton contribution from each of the lower latitude 
regions. In each frame of Fig.[4| the dashed line denotes 
this predicted contribution from inverse Compton scat- 
tering. Not surprisingly, this component makes up most 
or all of the total observed emission at high latitudes, 
but only a small fraction at low latitudes. Note that the 



^ The conclusion that no spectrum of electrons can produce the 
low-latitude gamma-ray emission could, in principle, be evaded 
if we had adopted a much more finely binned electron spectrum 
in our analysis. In particular, a cosmic ray electron spectrum 
that is described by a delta function at 16 GeV provides a good 
fit to the low-latitude feature. As such a feature is not expected 
from the perspective of cosmic ray acceleration, nor realistic in 
light of non- negligible energy loss processes, we do not consider 
this possibility further. 
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FIG. 7: The gamma-ray spectrum of the Fermi Bubbles after subtracting a contribution from inverse Compton emission, 
derived using the electron spectrum (up to normalization) found in our best-fit to the |6| = 40° — 50° region. This illustrates 
the characteristics of the additional (non-inverse Compton) component of the gamma-ray emission from the Fermi Bubbles, 
which is quite bright at low Galactic latitudes. We caution that these extracted spectra are subject to a number of systematic 
uncertainties, such as those associated with the interstellar radiation field model, and due to uncertainties and variations in the 
electron spectra throughout the volume of the Bubbles. These extracted spectra can, however, be taken as indicative of the 
broad spectral features of the non-inverse Compton component of the Bubbles emission. Shown as dashed lines is the predicted 
contribution of gamma-rays from the annihilations of 10 GeV dark matter particles (to r^r~) distributed according to a 
generalized NFW profile with an inner slope of 7 = 1.2, as described in Sec.jV] We remind the reader that the backgrounds are 
largest near the disk and thus there are significant systematic uncertainties in the spectrum from the low latitude (|6| = 1° — 10°) 
region, especially at low energies. 
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shapes of these dashed Hnes differ with latitude only as 
a result of variations in the radiation field model. 

To better isolate the spectrum of the non-inverse 
Compton emission, we subtract this estimate of the in- 
verse Compton contribution from the total emission of 
the Fermi Bubbles in each latitude range. This residual 
spectrum is plotted in Fig. [7] While one should keep in 
mind that the error bars shown in this figure do not take 
into account any variations in the spectral shape of cos- 
mic ray electrons throughout the volume of the Bubbles, 
this provides us with what is likely to be a reasonable es- 
timate of the spectrum and intensity of the non-inverse 
Compton emission exhibited in the low-latitude regions 
of the Fermi Bubbles. 

The spectra shown in Fig. [7] exhibit some rather dis- 
tinctive features. In particular, this spectral component 
peaks strongly at energies of ~l-5 GeV, and has no statis- 
tically significant presence above ~10 GeV. Furthermore, 
the intensity of this component is a very strong function 
of Galactic latitude, being more than an order of magni- 
tude brighter at low latitudes than at intermediate lati- 
tudes. These spectral and morphological characteristics 
are quite similar to those exhibited by the gamma-ray 
emission previously observed from the inner few degrees 
surrounding the Galactic Center [T5H2^ . In the following 
section, we will discuss this comparison in more detail. 



V. COMPARISON WITH GAMMA-RAY 
EMISSION FROM THE GALACTIC CENTER 

In previous studies of Fermi data, multiple sets of au- 
thors have identified the presence of a bright and spa- 
tially extended gamma-ray source around the Galactic 
Center, peaking at energies of a few GeV [THHH]. In 
particular, Ref. [HI reports that the morphology of this 
source implies a luminosity per volume that scales as 
r~'^-'^ to r~^ ®, where r is the distance to the Galactic 
Center. Similar profiles were found to provide good fits 
in Refs. [T51[?I1[^ . Each of these studies also found that 
the spectrum of this spatially extended emission peaks 
strongly at energies of a few GeV, very similar to the peak 
found in this study in the low-latitude emission from the 
Fermi Bubbles. 

This comparison strongly suggests that the non-inverse 
Compton emission we have identified in the low-latitude 
regions of the Fermi Bubbles is in fact the more spa- 
tially extended counterpart of the gamma-ray signal pre- 
viously reported from the innermost few degrees around 
the Galactic Center. To further explore this comparison, 
we plot as a dashed line in each frame of Fig.[7]the contri- 
bution predicted by an annihilating dark matter model 
found in Ref. [32, to provide a good fit to the Galactic 
Center signal. In this model, 10 GeV dark matter parti- 
cles annihilate to tau lepton pairs (with av = 2 x 10"^'' 
cm'^/s to T~^T~, OT av — 6 X 10"^'' cm^/s if annihila- 
tions proceed equally to r+T~, and e+e"), and are 
distributed following a generalized Navarro-Frenk- White 



profile, with an inner slope of 7 =1.2 (the classic NFW 
profile corresponds to an inner slope of 7 = 1), a scale 
radius of 20 kiloparsecs, and normalized such that the 
local dark matter density is 0.4 GeV/cm"^. Specifically, 
the functional form is: 

p(r) cx , 7 = 1.2, i?, = 20 kpc. (1) 

From this comparison of both the spectrum and mor- 
phology of these signals, we conclude that the gamma- 
ray component we have identified within the low-latitude 
regions of the Fermi Bubbles is the more spatially ex- 
tended continuation of the gamma-ray signal previously 
observed in the Galactic Center. Furthermore, we con- 
firm that the dark matter models previously shown to be 
capable of accounting for the gamma-ray emission in the 
inner Galaxy are also capable of providing an explana- 
tion for the non-inverse Compton emission observed in 
the low-latitude spectrum of the Fermi Bubbles. 

VI. MORPHOLOGY OF THE LOW-LATITUDE 
SIGNAL 

One approach to extracting the morphology of the few- 
GeV spectral feature discussed in the previous sections 
is to examine the residual sky maps produced when the 
background templates (for example, the diffuse model 
template and isotropic background), multiplied by their 
best-fit coefficients, are subtracted from the data. Equiv- 
alently, these maps are produced by taking the residu- 
als of the template fit, and re-adding the latitude-sliced 
Bubbles templates, weighted by their best-fit coefficients. 
This largely cancels out structure due to mismatches be- 
tween the shape of the Bubbles templates and the actual 
excess, although such mismatches may still bias the fit. 
(This procedure has been recommended for extracting 
the spectrum of an excess in Ref. |/7j.) 

Fig. |8] shows these "residual" maps (in average 
E^dN/dE) for the "diffuse model" fit, masked within 5 
degrees of the Galactic plane, in four energy bands span- 
ning the range from 1 — 20 GeV. The high-latitude Fermi 
Bubbles are visible in all bins with comparable bright- 
ness, as expected since the high-latitude spectrum is ap- 
proximately flat in E'^dN/dE. At low latitudes, there 
is a pronounced excess around the Galactic Center, not 
disk-like in shape, which is clearly visible in the 1 — 2 and 
2 — 5 GeV maps, but rapidly fades at higher energies. 

Having generated these residual maps, we can study 
the emission binned in &, within a narrow longitude range 
(1^1 < 5). As an illustrative exercise, we subtract the 
10 — 50 GeV (residual) map from the 1 — 10 GeV (resid- 
ual) map, where the excess is most pronounced (we ignore 
energies higher than 50 GeV because the P6V11 Fermi 
diffuse model was not designed for studies at those en- 
ergies and, due to low photon statistics, little informa- 
tion would be gained from higher-energy photons in any 
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FIG. 8: The residual emission after re-adding the latitude-shced Bubbles templates with their best-fit coefficients, in E^dN/dE. 
Equivalently, these maps are obtained by subtracting the best-fit model for the background (in which we include all templates 
but the latitude-sliced Bubbles) from the data. The "diffuse model" fit is used, performed over regions greater than 5° from 
the plane (although the mask shown in the figure is at |6| =3°). 
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FIG. 9: The difference of the residual emission maps (after re-adding the latitude-sliced Bubbles templates with their best-fit 
coefficients), between the 1-10 GeV bin and the 10-50 GeV bin, in E^dN/dE averaged over —5° < I < 5° . The error bars 
describe the pixel-to-pixel scatter within each bin (standard deviation of pixel values). This analysis employs the "diffuse 
model" fit (see text), masked at 5° from the plane. The red line shows the anticipated intensity resulting from a (squared, 
projected) NFW profile with inner slope of 7 = 1.2. 



case). For more details, see Appendix |e] Since the maps 
are given in E'^dN/dE, a zero result indicates an average 
spectrum of dN/dE oc E~'^ between these two bins. The 
results for the southern hemisphere, where there are fewer 
bright sources and local features, are shown in Fig. |9] 
While we find that at high latitudes the spectrum is con- 
sistent with dN/dE oc E~'^, the lower latitude emission 



{\h\ < 10° — 15°) reveals significant additional emission 
at low energies. The error bars shown in this figure, com- 
puted from the standard deviation of the pixel values in 
each Ab = 2° bin, are quite large and non-negligibly cor- 
related, but provide a sense of the uncertainty in the rate 
at which the signal falls off away from the Galactic plane. 
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FIG. 10: The spectrum of the emission fit by the NFW template (7 — 1.2). The far left, center left and center right panels 
employ the "diffuse model" fit, while the far right panel uses the "low-energy template" fit (see text for details). In all cases 
an additional (squared, projected) NFW template is added to the fit, and its extracted spectrum is plotted. The center left 
panel also includes a template for Loop 1, and in the center right panel the fit is performed using the southern hemisphere 
only, where there is less bright residual structure. The different colors show different choices for the latitude cut to remove the 
Galactic Disk: |b| < 1° (black), |6| < 2° (blue), |ti| < 5° (red). Where the Icr error bars overlap with zero, we instead plot 
downward-pointing arrows corresponding to the 3(t upper limits on the emission. 



VII. FITTING THE SPECTRAL BUMP WITH 
AN NFW PROFILE 

So far, we have focused our attention on the region of 
the sky occupied by the Bubbles. If annihilating dark 
matter is responsible for the signal observed at low lati- 
tudes, however, this signal should be distributed around 
the Galactic Center with approximate spherically sym- 
metry, and will not be restricted to the Bubbles. In this 
section, we include an additional template in our anal- 
ysis, corresponding to the signal expected from annihi- 
lating dark matter distributed according to a somewhat 
steepened NFW profile (with an inner slope of 7 = 1.2). 

In Figs. [T0| and [Tl] we show the spectra associated with 
the NFW template itself, and with the latitude-sliced 
Bubbles in the presence of the NFW template, respec- 
tively. The normalization of the NFW template corre- 
sponds to its brightness at 5° from the Galactic Center. 

With the inclusion of this additional NFW template, 
we find that the low-latitude (non-inverse Compton) 
emission prefers the morphology of the NFW template 
over being absorbed into the low-latitude portion of the 
Bubbles. The significant improvement in the fit is driven 
both by the spherical symmetry of the NFW profile, and 
the fact that the imposed Bubbles templates are fiat in 
brightness rather than falling off with increasing Galacto- 
centric radius (see also Appendix [d|) . With the inclusion 
of this NFW template in the fit, the spectrum of the 
emission associated with the Bubbles is not only largely 
constant at high latitudes (as expected), but also has 
a similar (near-flat) spectral shape at low latitudes, at 
least where it can be detected (in several cases there is 
no significant detection of Bubbles-correlated emission 
in the lowest-latitude band, after the addition of the 
NFW template). The significant GeV-scale feature iden- 
tified in Sec. |IV| is no longer absorbed by these Bubbles 
templates, but instead is present in the spectrum of the 
spherically symmetric NFW template. The spectrum as- 



sociated with the NFW template peaks at ~l-2 GeV, 
and falls off above ~10 GeV (there is little or no evi- 
dence for emission above 10 GeV in the spectrum of the 
NFW template, at least in the diffuse model fit). The 
total flux (integrated over solid angle) associated with 
this template is ^ 1 keV/cm^/s, corresponding to a total 
luminosity within the solar circle of ~ 10^^ erg/s. 

The regions of the sky that are most important for dis- 
criminating between the NFW template and the Bubbles 
templates are those near the inner Galaxy but outside the 
Bubbles themselves. Unfortunately, such regions are also 
near the Galactic plane. As a result, the extracted spec- 
trum associated with the NFW template depends some- 
what on the details of the fit, including the mask of the 
Galactic plane, especially below ~1 GeV where the an- 
gular resolution of Fermi is somewhat poor. In the sim- 
plest case, in which we employ the "diffuse model" fit and 
add the NFW template, the spectrum we extract below 
1 GeV is considerably softer than that extracted from 
the Bubbles templates (see the left frame of Fig. 10 1. 



Masking |&| < 5° hardens the low-energy spectrum, as 
does adding a separate template to absorb soft-spectrum 
emission from Loop I, or restricting the fit to the southern 
hemisphere where the structures associated with Loop I 
do not contaminate the signal (see left center and right 
center frames of the figure). Given these results, we find 
it plausible that the low-energy emission (below 1 GeV) 
associated with the NFW template can be attributed in 
large part to contamination from the disk and Loop I. 

In Fig. [12] we show how the spectral feature at a few 
GeV is reallocated from the lowest-latitude slices of the 
Bubbles to the NFW template, once the NFW template 
is added. In order to minimize background contamina- 
tion, we examine the spectra masking |6| < 5° and re- 
stricted to the southern hemisphere. As discussed above, 
we find that the "bump" is much better fitted by the 
NFW template and there is no remaining significant evi- 
dence for emission correlated with the low-latitude Bub- 
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FIG. 11: The spectrum extracted for the gamma-ray Bubbles in ten-degree latitude bands, now including an NFW template 
in the fit: in order from the top row, 40° < |6| < 50°, 30° < |b| < 40°, 20° < |fe| < 30°, 10° < |6| < 20°, and |6| < 10°. The left 
and center panels use the "diffuse model" template fit (see text); in the center panels, the fit is restricted to & < in addition 
to the masking. The right panels use the "low-energy template" fit (see text). In all cases an additional (squared, projected) 
NFW template is added to the fit, as in Fig. |10[ The different colors show different choices for the latitude cut to remove the 
Galactic Disk: |6| < 1° (black), |fe| < 2° (blue), |6| < 5° (red). Where the Icr error bars overlap with zero, we instead plot 
downward-pointing arrows corresponding to the 3a" upper limits on the emission. 
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FIG. 12: The spectra of the various fit components, including five separate latitude-sliced templates for the Bubbles (see text), 
for the diffuse model foreground template, when the fit is restricted to the southern sky only. The Galactic Disk is masked 
for |6| < 5° in all cases. The left and right panels correspond to the case with (right) and without (left) an additional NFW 
(7 — 1.2) template in the fit. 
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FIG. 13: The formal statistical -~2Aln£, values (referred to as A^^) extracted from the likelihood fit using the diffuse model, 
the sliced Bubbles templates and a projected squared NFW template, for varying values of the NFW inner slope 7. The Ax^ 
values are with respect to the tested 7 value with the highest likelihood (so the minimum is zero by definition). Red stars 
correspond to the full-sky fit masked for j6j < 1°, black diamonds to the southern-sky fit masked for |6j < 5°. 
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bles. Quantitatively, we find that the fit including the 
NFW template improves over the five-Bubbles-templates 
model at the level of approximately 12a. Again, though, 
we reiterate that this is a formal significance derived in- 
cluding statistical errors only, whereas there are signifi- 
cant systematics originating from the fact that the mod- 
els we are using do not fully describe the data. 

We have performed fits using generalized NFW pro- 
files with a variety of inner slopes, 7. Fitting over both 
hemispheres and masking within 1 degree of the Galactic 
plane, we find that 7 ~ 1.2 provides the best fit (for this 
reason, we have adopted 7 =1.2 as our default value). If 
we restrict our fit to the southern hemisphere and mask 
within five degrees of the Galactic plane, however, we 
find that a somewhat steeper distribution is preferred, 
7 ~ 1.5 — 2. This may refiect a slope which varies with 
distance from the GC, although systematic uncertainties 
prevent us from making any strong statement to that ef- 
fect. We show the Aln£ values for various choices of 7, 
in both these cases, in Fig. 
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We note that it is not surprising that a spherically 
symmetric signal centered around the Galactic Center 
would first become apparent in an analysis of the region 
of the Fermi Bubbles. By their shape, the Bubbles ex- 
clude both most of the background from the disk, and 
the arc structures associated with Loop I, making them 
reasonably well designed for extracting even a spherically 
symmetric excess (note the similarities between the Bub- 
bles regions and the dark matter regions-of-interest as 
described in Refs. 



VIII. INTERPRETATION 

In the existing literature, a number of possibilities have 
been discussed for the origin of the extended gamma-ray 
emission observed from the region of the Galactic Center, 
including annihilating dark matter [T8'-'20l [22] , a popula- 
tion of millisecond pulsars [HI [19, 22 24j , and cosmic ray 
interactions with gas [HI [IH [HI 1251 . In this section, 
we discuss these possibilities in light of the new informa- 
tion provided in this study, including the evidence that 
this emission is not confined to the inner Galaxy, but in- 
stead extends out to at least 2-3 kpc from the Galactic 
plane. 



Diffuse Emission Mechanisms 



As we have previously argued in Sec. |IV[ no realistic 
spectrum of cosmic ray electrons can produce the ob- 
served spectral features of the GeV emission identified 
in this study. Cosmic ray protons scattering with gas 
also fail to provide a reasonable explanation for the ob- 
served gamma-rays. Firstly, the morphology of the ob- 
served emission is not highly correlated with the distri- 
bution of gas in the Milky Way, as would be predicted 
in such a scenario. Secondly, the spectral shape of the 



observed gamma-ray emission requires a peculiar spec- 
trum of cosmic ray protons, sharply peaked at approx- 
imately 25 GeV. In the upper left frame of Fig. [T4j we 
compare the observed spectrum of this emission to the 
spectral shape predicted from proton collisions with gas 
for a proton spectrum given by a delta function at 25 GeV 
(solid), or by a broken power-law, following E~'^ below 
25 GeV and E^^ at higher energies (dashed). While nei- 
ther provides a particularly good fit, it is clear that from 
this comparison that a very sharply peaked spectrum of 
protons would be required to potentially account for this 
emission. And while we cannot absolutely rule out the 
existence of such a strongly peaked cosmic ray proton 
spectrum, this rather extreme requirement further dis- 
favors cosmic ray origins for the observed low-latitude 
emission. 



Note that in each frame of Fig. 14 we make compar- 
isons to the observed gamma-ray spectrum (after sub- 
tracting the contribution from inverse Compton), as ex- 
tracted using the diffuse model template and from the 
\b\ = 10° - 20° regions of the Bubbles. We have cho- 
sen to compare to this region over that extracted from 
lower latitudes, or from the NFW-template, as this re- 
gion provides a spectrum that is more robust to contam- 
ination from the Galactic Disk. In particular, a com- 
parison of the low-energy spectrum extracted from the 
|6| < 10° Bubbles template or (especially) the NFW tem- 
plates varies considerably depending on how the disk is 
masked, and on whether we consider both hemispheres, 
or only the southern sky. On the other hand, the spec- 
tral shape extracted at higher latitudes depends more 
strongly on the subtraction of the inverse Compton com- 
ponent, which is negligible compared to the bump for 
\b\ < 10°, but increasingly substantial as one moves far- 
ther from the Galactic plane. However, the consistency in 
the spectrum of the bump between the 10—20° band and 
even higher latitudes supports the simple picture where 
the electron spectrum generating the ICS component (or 
the proton spectrum, in a hadronic scenario) is essentially 
invariant. 



B. Millisecond Pulsars 

Instead of relying on diffuse emission mechanisms, a 
large population of unresolved gamma-ray point sources 
could, in principle, be responsible for the observed emis- 
sion. In particular, the spectrum of gamma-ray pulsars 
has been observed to peak at GeV energies, and it has 
been suggested that a collection of ~10"^ such objects may 
be able to account for the extended gamma-ray emission 
observed from the Galactic Center [22l - f24] . Of partic- 
ular interest are millisecond pulsars, which are thought 
to be formed as parts of low-mass X-ray binary systems. 
This could help to accommodate the very concentrated 
distribution of the gamma-ray emission observed around 
the Galactic Center (the stellar distribution in the inner 
Galaxy scales as n c>c r~^-'^, but objects formed through 
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FIG. 14: Comparisons of the observed gamma-ray spectrum of the low-latitude emission, after subtracting the contribution 
from inverse Compton scattering (see Fig.[7]l to that predicted from the scattering of cosmic ray protons with gas {upper left), 
millisecond pulsars (upper right), and dark matter annihilations (lower). For proton-gas collisions, the solid and dashed lines 
denote cosmic ray proton spectra which take the form of a delta function at 25 GeV or a broken power-law following Ep^ below 
25 GeV and Ep'^ at higher energies, respectively. To accommodate the spectral shape of the observed gamma-ray emission, 
the cosmic ray proton spectrum throughout the inner several kiloparces of the Fermi Bubbles must peak very strongly at 
approximately 25 GeV. The spectrum shown for pulsars is that corresponding to the average millisecond pulsar observed by 
the Fermi collaboration |36l 137] . For annihilating dark matter, we show results for two models: 10 GeV particles annihilating 
to tan leptons (dashed) and 50 GeV particles annihilating to bb. In each case, we have adopted a generalized NFW profile with 
an inner slope of 7 = 1.2, and normalized the signal to a local density of 0.4 GeV/cm'^ and an annihilation cross section of 
av = 2x 10"^'^ cm^/s (r+r~) or cru = 8 x 10"^^ cm^/s (bb). 



interactions between stars could plausibly be distributed 
as steeply as the square of this distribution) [22l|38]. The 
binary companions of such pulsars could also act as a 
tether, explaining why they do not free-stream out of the 
Galactic Center as a result of velocity kicks. Further- 
more, millisecond pulsars are spun up through accretion, 
and can thus produce high luminosities of gamma-ray 
emission over much longer timescales than other types of 
pulsars. 

Unfortunately, relatively little is known about the 
gamma-ray spectrum of millisecond pulsars. Fermi has 
reported spectra from only eight millisecond pulsars, 
which together yield an average spectrum that is well fit 
by dN^/dE^ cx E-'^-^exp{-E^/2.8GeV) giEZ]. In the 



right frame of Fig. 14 we show this spectral shape (with 



arbitrary normalization) compared to the observed emis- 
sion from the low-latitude Bubbles. This does not pro- 
vide a good fit, especially at low energies (although this is 
also where the potential systematic errors are most signif- 
icant). One could imagine, however, that the eight mil- 
lisecond pulsars reported by Fermi may not be represen- 
tative, perhaps being biased toward the brightest or most 
locally common examples of such objects. To address this 
issue, it has been suggested that the gamma-ray spectra 
of globular clusters (which are thought to contain large 
numbers of millisecond pulsars) could provide a more re- 
liable determination of the average spectrum from mil- 
lisecond pulsars 23 . At present, however, the error bars 
on the gamma-ray spectra of globular clusters are quite 
large, and (on average) do not appear to favor a much 
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harder spectral index than is observed from individual 
resolved pulsars [39] (see also Fig. 9 of Ref. |l9j). Future 
data from Fermi could potentially be useful in further 
testing the possibility that millisecond pulsars produce a 
gamma-ray spectrum compatible with the low-latitude, 
GeV emission under consideration here. 

In addition to their gamma-ray spectra, the distribu- 
tion of millisecond pulsars within the Milky Way is not 
well constrained empirically. That being said, one ex- 
pects the formation of millisecond pulsars to roughly fol- 
low that of stars, possibly with an additional preference 
for regions of very high stellar density, such as in globu- 
lar clusters. To account for the observed morphology of 
this signal, however, there would have to be a significant 
number of millisecond pulsars in the halo, at distances 
of at least a few kpc outside the Galactic plane. If this 
emission is generated by such objects, they would require 
a distribution that is very different from that observed 
among other stellar populations. And while such distri- 
butions have been proposed [40], such distributions are 
constrained by the small number of millisecond pulsars 
resolved by Fermi and by the high degree of isotropy 
observed in the gamma-ray background |41j . 



C. Annihilating Dark Matter 

In the previous two sub-sections, we have described 
some of the possibilities and challenges involved in ex- 
plaining this gamma-ray signal with astrophysical sources 
or mechanisms. Annihilating dark matter can provide 
a simple explanation for the sharply peaked spectrum 
and distinctive morphology of this emission. In the lower 
frame of Fig. |14[ we compare the gamma-ray spectrum 
predicted in two dark matter scenarios to the observed 
spectrum. First, we consider a 10 GeV dark matter par- 
ticle species annihilating to tau leptons. We also show 
the spectrum resulting from a 50 GeV particle annihilat- 
ing to bb which, given the systematic uncertainties in the 
extraction of the spectrum, we also consider to be a vi- 
able possibility. For a generalized NFW profile (7 = 1.2) 
normalized to a local density of 0.4 GeV/cm'^, we require 
an annihilation cross section of av = 2 x lO"^*" cm^/s 
in the 10 GeV t+t^ case, or av — 8 x 10"^'' cm^/s 
in the case of a 50 GeV particle annihilating to bb. As 
annihilations to electrons, muons, or neutrinos do not 
significantly contribute to the gamma-ray spectrum, the 
total annihilation cross section in the leptonic case may 
be larger by a factor of a few or several. 

We note that the dark matter distribution required to 
fit the observed signal is well supported by the results of 
numerical simulations. As has been known for some time, 
N-body simulations which model the evolution of cold 
dark matter without baryons typically find halos of the 
NFW form, with inner slopes oi p (x (7 = 1) [42ll43j . 
Hydrodynamical simulations, which include the effects of 
baryonic processes involved in galaxy formation and evo- 
lution, in many cases predict the steepening of this inner 



slope, typically from 7 = 1 to 7 = 1.2 to 1.5 (see Ref. [ii] 
and references therein). Supporting this picture further 
are observations of the Milky Way's rotation curve and 
microlensing optical depth, which together find that the 
dark matter distribution is best fit by 7 = 1.3 (although 
with large uncertainties) ^45j.'* 

The annihilation cross section required to normalize 
the signal in question is also an attractive value from 
a theoretical perspective. To be produced thermally in 
the early universe in an abundance equal to the mea- 
sured dark matter density, a (~10-50 GeV) dark mat- 
ter particle must have an annihilation cross section of 
{(jv) K, 2 X 10~^^ cm^/s, when thermally averaged over 
the process of freeze-out [48 . In many typical dark mat- 
ter models, however, the annihilation cross section in the 
low velocity limit (as is relevant for annihilation in the 
halo) is smaller than this value by a factor on the or- 
der of a few. Neutralinos, for example, when annihilat- 
ing to leptons through the i-channel exchange of slcp- 
tons typically exhibit (crw)Fo/o"^ii=o ^ 5. While neu- 
tralino annihilations through p-wave amplitudes (which 
are suppressed at low velocities) significantly contribute 
to the relic abundance calculation, they do not factor 
into the current annihilation rate. We also note that 
annihilation cross sections in the range required here 
[av ~ (2 — 8) X 10^^^ cm'^/s) are currently consistent 
with all observational and experimental constraints, in- 
cluding those derived from dwarf galaxies [49l [50], the 
cosmic microwave background [51H53j . the spectrum of 
cosmic-ray antiprotons [541 I55j . and from searches for 
monophoton-plus-missing energy events at LEP [5^ . 

From a dark matter model building perspective, there 
are a number of possibilities one could consider. In 
the case of light dark matter particles which annihi- 
late primarily to leptons, the annihilations could proceed 
through the t-channel exchange of new charged parti- 
cles which carry lepton number (such as sleptons), or 
thorough the exchange of a leptophillic Z' [57] or a lep- 
tophillic higgs [58]. Alternatively, one could imagine a 
scenario in which the dark matter itself carries lepton 
number. The dark matter could also annihilate to a pair 
of light force carriers which interact with the Standard 
Model only through kinetic mixing the with photon. The 
decays of the light force carrier then yield combinations 
of mesons and charged leptons, leading to a gamma-ray 
spectrum that is very similar to that predicted from tau 
decays [59] . 

One additional consideration is that if the ~10 GeV 
dark matter particles also annihilate to e^e~ at a rate 
similar to r+r^, then they would also be capable of ex- 
plaining the anomalous synchrotron emission observed 
from the Milky Way's radio filaments [50], and possi- 



* Note, however, that the inclusion of baryonic feedback processes 
may flatten cusps in dark matter halos (e.g. 1461 ). and there is 
observational evidence for profiles shallower than NFW in some 
dwarf galaxies |47| . 



19 



bly a significant fraction of the isotropic radio back- 
ground [611 E2] ■ Such annihilations could also lead to a 
feature in the cosmic ray positron spectrum, potentially 
identifiable by AMS [53]. 

IX. DISCUSSION AND CONCLUSIONS 

In this study, we have identified a new component of 
gamma-ray emission from the low-latitude regions of the 
Fermi Bubbles. This emission does not appear to be 
compatible with originating from either inverse Compton 
scattering or proton collisions with gas (for any physi- 
cally realistic spectrum of cosmic rays). The spectrum 
of this emission peaks at energies of a few GeV and is 
distributed with an approximately spherically symmetric 
morphology about the Galactic Center, with a luminosity 
per volume that is consistent with a squared (generalized) 
NFW profile with an inner slope of p (x r~^-^ (although 
the best-fit slope depends somewhat on the region chosen 
for the fit). The broad features of this signal are robust 
to the various details of our analysis, such as the choice 
of templates and the degree of masking around the disk. 

The diffuse gamma-ray signal we have identified here is 
consistent both in spectrum and morphology with being 
the more spatially extended counterpart of the emission 
previously observed from the Galactic Center [T5ti22j . 
In addition to further confirming the existence of this 
gamma-ray source, the results presented here are impor- 
tant in that they confirm that the morphology of this 
signal extends to at least ^2 — 3 kpc from the Galactic 
Center (whereas previous studies of the Galactic Center 
were sensitive only to the inner few hundreds of parsecs) . 
This provides us with valuable information with which 
to test various interpretations of the observed emission. 
Millisecond pulsars, in particular, have been proposed 
as a source of the emission observed from the Galactic 
Center. In such a scenario, the significant high-latitude 
component of the signal reported in this paper would 



seem to require a large high-latitude population of faint 
unresolved millisecond pulsars that does not trace known 
stellar populations. In contrast, the morphology of the 
signal is consistent with that predicted a priori from an- 
nihilating dark matter. 

If annihilating dark matter is responsible for the emis- 
sion reported here, we find that a halo profile of the gen- 
eralized NFW form, with an inner slope of p oc r"^'^, is 
preferred. The spectral shape of the emission is well de- 
scribed by an approximately 10 GeV particle annihilating 
to tau leptons, with a normalization corresponding to an 
annihilation cross section of cru ^ 2 x 10~^^ cm'^/s (in 
addition to annihilations proceeding to electrons, muons 
and/or neutrinos, each of which provide comparatively 
few gamma-rays). We note that dark matter in this mass 
range may also be capable of accommodating a number 
of other anomalous signals reported from direct and indi- 
rect detection experiments (for a summary, see Ref. 32]). 
Alternatively, the annihilations of a ^-^40-50 GeV dark 
matter particle to quarks can also provide a reasonable 
fit, requiring a normalization of crv ^ 8 x 10^^*" eve? /s. 
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Appendix A: The Simple Disk model 

The primary template subtraction technique initially 
used in Ref. \Xi to reveal the Bubbles employed a sim- 
ple geometric template to model ICS emission from the 
disk. The functional form of this template is esc \b\ — 1 in 
latitude and Gaussian with cr; — 30° in longitude. Such 
a model, while oversimplified relative to the actual emis- 
sion associated with the Galactic Disk, is smooth and 
cannot mimic sharp features, and is demonstrably effec- 
tive in subtracting ICS emission at high latitudes. 

In the current work, we attempted to fit the data us- 
ing this template to model ICS, along with the SFD 
dust map to trace emission from pion decay, the isotropic 
background model, a flat template for Loop 1, and the 



latitude-sliced Bubbles templates. This is identical to 
the "low-energy template" fit described in the main text 
except that the dust-subtracted 0.5 — 1 GeV map is re- 
placed with the simple disk. The results of this approach 
are shown in Fig. |15| This model turns out to be poorly 
suited to fits at low latitudes, with the fit being highly 
unstable to the masking of the disk; we suspect that 
the steep rise of the template at small \b\ causes the fit 
to be largely driven by the lowest-latitude data, where 
the data are known to be contaminated by unsubtracted 
point sources. 
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FIG. 15: The spectrum extracted for the gamma-ray Bubbles in ten-degree latitude bands; in order from the top left, 
40° < \b\ < 50°, 30° < \b\ < 40°, 20° < |&| < 30°, 10° < Ife] < 20°, < 10°. The fit is performed using the simple disk 
model as a template for ICS, the SFD dust map, a uniform map, a template for Loop I, and the five latitudinal slices through 
the Bubbles templates. The different colors show different choices for the latitude cut to remove the Galactic Disk: |6| < 1° 
(black), l&l < 2° (blue), |6| < 5° (red). Where the la error bars overlap with zero, we instead plot downward-pointing arrows 
corresponding to the 3a upper limits on the emission. The bottom right panel shows the spectrum of the various templates in 
the fit for the |6| < 1° mask. 



Appendix B: Consistency of results using the diffuse 
model and low energy template 

At first appearance, the gamma-ray spectra as ex- 
tracted using the diffuse model and low-energy template 
methods are quite different (see the left and right frames 
of Fig. [3|. This comparison can be misleading, however, 
as the spectra extracted using the low-energy template 
do not describe the emission present in the Bubbles re- 
gion after subtraction of some physical model - rather, 
they describe the degree to which the Bubbles-correlated 
emission is harder than the ordinary ICS emission as- 
sociated with the disk (which dominates the 0.5-1 GeV 
skymap after removal of the dust-correlated emission), 
with a pivot point at 0.5-1 GeV. In contrast, the diffuse 
model template fit extracts the excess over the diffuse 
model spatially correlated with the Fermi Bubbles. Of 
course, the diffuse model is fitted to the data, and so 
may also soak up emission physically associated with the 



Bubbles. 

In this Appendix, we provide a direct comparison of 
these techniques. Comparing the gamma-ray spectra ex- 
tracted using the two methods, we see that there are dis- 
crepancies at both low and high energy (with the latter 
being particularly pronounced at low latitudes). These 
discrepancies have different origins; we will discuss them 
in turn. 

At high energies and low latitudes, the low-energy tem- 
plate fit returns a roughly flat (in E^dN/dE) Bubbles- 
correlated spectrum, similar to the spectrum obtained 
at high latitudes. In contrast, the diffuse model fit de- 
tects virtually no emission. This indicates that the dif- 
fuse model is "soaking up" this hard-spectrum emission 
in some way. The diffuse model has, of course, been 
adjusted to fit the data, and has the potential to over- 
subtract Bubbles-correlated emission in trying to fit the 
data without the freedom to include the Bubbles explic- 
itly. Given the information to hand, it is not possible 
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to distinguish between this scenario and one where the 
Fermi diffuse model is adequately capturing the physics 
of this spectral component, which should be assigned to 
the ordinary Galactic emission even though it is harder 
than the norm. The amplitude of this hard spectral com- 
ponent, at low latitudes, is not of great interest for the 
purposes of this study; to facilitate checking agreement of 
the few-GeV spectral feature between the two methods, 
we subtract a component with dN^/dE^ cx E^"^ from 
the results for the low-energy template fit, designed to 
remove this high-energy emission. 

The low-energy discrepancy is simpler to deal with, 
as it is an entirely natural and expected consequence of 
the different methodologies being employed. Since the 
low-energy template already includes some contribution 
from the Bubbles, to obtain the true Bubbles spectrum 
we should re-add a component with the same spectrum as 
the low-energy template {dN^/dEj cx E~'^-^^, to a good 
approximation). The normalization of this component is 
a priori unknown, but to check consistency between the 
two methods, we can normalize it to match the diffuse 
model result at 0.5-1 GeV. 

Having chosen these two free parameters (the ampli- 
tude of the high-energy spectrum and the amplitude of 
the correction due to the low-energy template), we com- 
pare in Fig. [16] the spectra extracted using the diffuse 
model template (black) to that extracted using the low- 
energy template, after performing the adjustments de- 
scribed in this Appendix (blue). After these corrections, 
the low-energy template results are almost identical to 
those derived using the diffuse model template. 



Appendix C: The symmetry of the signal 

One might ask if the signal possesses north-south sym- 
metry, as well as left-right symmetry (negative vs positive 
I). In Fig. 17 we show the variation in the results when 



the fit is performed only in the north, only in the south, 
only for ? < 0, or only for I > 0. In all cases, the spectral 
feature is apparent in the lowest latitude bin. The largest 
difference is in the relatively flat-spectrum (in E'^dN/dE) 
emission attributed to the high-latitude Bubbles; this is 
larger in the north and for I > 0. This is to be expected, 
as Loop I and the associated gamma-ray arcs occur in 
the north quadrant with I > 0. 



Appendix D: Searching outside the Bubbles 



In Sec. |VII[ we showed that the low-latitude GeV-scale 
spectral feature is considerably better fit by a spherically 
symmetric NFW profile than by the Bubbles templates. 
In this Appendix, we perform a cross check of this conclu- 
sion by explicitly asking whether this spectral component 
is observed from the regions of the sky not associated 
with the Bubbles. The absence of such emission from 



the "bubble complement" region would strongly disfavor 
any dark matter interpretation of the signal. 

To address this question, we re-fit the data with an 
additional template, corresponding to the region lying 
within 20° of the Galactic Center, but not within the 
Bubbles. This region is indicated by dotted lines in 
Fig. [T] Fig. [18] shows the results for the "diffuse model" 
and "low-energy template" fits with this component 
added. 

In each frame (corresponding to differing degrees of 
masking of the disk), the spectrum of the Bubbles 
complement template exceeds or equals the predicted 
gamma-rays from dark matter (for the same dark mat- 
ter model and distribution as shown by dashed lines in 
Fig. [t]). In contrast, little or no such emission is associ- 
ated with this template at energies above 10 GeV. 

At low energies, the spectrum and amplitude of the 
Bubbles complement template is quite dependent on the 
degree to which we mask the Galactic Disk, and has re- 
semblances in shape to the spectra associated with the 
NFW template. This supports our earlier hypothesis that 
the low-energy emission associated with the NFW tem- 
plate in these cases is being largely driven by the struc- 
tures associated with Loop I and the Galactic disk. 



Appendix E: Residuals 



As discussed in Sec. |VI| residual maps with the 
latitude-sliced Bubbles re-added provide a characteriza- 
tion of the excess emission which is less dependent on 
the choice of signal templates. As in the main text, we 
restrict ourselves to |^| < 5°, and bin in latitude in two 
degree steps, Ab — 2°. In Fig. [TOj we show the resulting 
"Bubbles residuals" for energy bands from 0.5-1 GeV, 1- 
10 GeV, and 10-50 GeV. The residuals are relatively flat 
in b for |&| > 10°, and similar (in E'^dN/dE) for the three 
energy bands. At smaller in contrast, the results for 
low energies diverge from those for the high-energy band. 

In Fig. [20j we look more closely at the 1 — 10 GcV 
energy band, where the signal is present. We display 
the total emission, the Bubbles residual (i.e. the residual 
of the fit -I- the best-fit Bubbles templates), and the fit 
residual. We see that, with the exception of the region 
close to the Galactic plane where the Bubbles templates 
vanish, the absolute value of the fit residual is generally a 
factor of a few lower than the best-fit Bubbles-correlated 
emission, and fluctuates between positive and negative 
values; the Bubbles-correlated emission is another factor 
of several below the total emission in this energy band, 
depending on latitude. Thus, except for the region clos- 
est to the Galactic plane, this approach appears to pro- 
vide an adequate model of the features in the data, and 
unaccounted-for residuals are unlikely to skew the signal. 
Close to the Galactic plane, there are large residuals, as 
the small extent of the Bubbles means they are unable 
to absorb an extended signal; given the uncertainty in 
the shape of the Bubbles close to the plane, this is not 
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FIG. 16: A comparison of the gamma-ray spectra extracted from various latitude regions of the Fermi Bubbles using the diffuse 
model template (black) to that extracted using the low-energy template, after performing the corrections described in the text 
(blue). After accounting for these corrections, the low-energy template results are almost identical to those derived using the 
diffuse model template. 



unexpected, and motivates the addition of a template to absorb this additional emission. 
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FIG. 17: Tlie spectra of the various fit components, including five separate latitude-sliced templates for the Bubbles (see text), 
for the diffuse model foreground template, restricting the fit to 6 < {top left), b > (top right), I < {bottom left), I > 
{bottom right). The Galactic Disk is masked for |6j < 1° in all cases. 
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FIG. 18: The spectrum of the "bubble complement" region, defined as y/b'^ + P < 20° and that does not lie within the 
Bubbles templates, as extracted from a fit including the complement template, the five latitudinally sliced component Bubbles 
templates, and background templates as described in Fig. [3] Here, we have used the Fermi diffuse model. The different frames 
correspond to different choices for the latitude cut to remove the Galactic disk: |6| < 1°, |6| < 2°, |fe| < 5°. (Note that this 
substantially changes the region over which the complement spectrum is averaged, and thus may truly modify the result.) The 
dashed line in each frame denotes the average spectrum expected in this region from the dark matter model and distribution 
as shown in Fig. [7] 
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FIG. 19: The residual emission after re-adding the latitude-shced Bubbles templates with their best-fit coefficients, as a 
function of energy. Emission is E^dN/dE averaged over \l\ < 5° and in bins of width A6 = 2°; the error bars describe the pixel- 
to-pixel scatter within each bin. The upper panel shows the fit using the diffuse model, the lower the fit using the low-energy 
template, both masked at 5° from the plane. 
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FIG. 20: The absolute value of the residual emission after subtraction of all templates (dotted line) and the emission after 
re-adding the latitude-sliced Bubbles templates with their best-fit coefficients (red stars), for 1-10 GeV. The total emission is 
shown in black diamonds. Emission is E^dN/dE averaged over |/| < 5° and in bins of width Ab — 2°; the error bars describe 
the pixel-to-pixel scatter within each bin. The second and fourth panels show the residual emission after subtraction of all 
templates (not its absolute value). The upper two panels show the fit using the diffuse model, the lower the fit using the 
low-energy template; in both cases the fit was performed for |6| > 5°. 



